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The crossover from nonadiabatic to adiabatic electron 
transfer has been theoretically studied under a spin-boson 
model (dissipative two-state system) description. We present 
numerically exact data for the thermal transfer rate and 
the time-dependent occupation probabilities in largely un- 
explored regions of parameter space, using real-time path- 
integral Monte Carlo simulations. The dynamical sign prob- 
lem is relieved by employing a variant of the recently proposed 
multilevel blocking algorithm. We identify the crossover 
regime between nonadiabatic and adiabatic electron trans- 
fer, both in the classical (high-temperature) and the quantum 
(low-temperature) limit. The electron transfer dynamics dis- 
plays rich behaviors, including multi-exponential decay and 
the breakdown of a rate description due to vibrational coher- 
ence. 



I. INTRODUCTION 

Electron transfer (ET) reactions are important for an 
understanding of many processes in chemistry, biology 
and physics. Popular examples include charge transfer in 
semiconductors, chemical reactions, or the primary ET 
step in bacterial photosynthesis, and, despite the long 
history of the field 0,0], ET still attracts unbroken and 
(probably exponentially) increasing attention In 
this paper, a theoretical study of ET reactions based on 
the spin-boson (dissipative two-state) description is pre- 
sented. The spin-boson model provides a well-established 
description of ET processes in condensed phase environ- 
ments, which can be motivated from the central limit the- 
orem or, alternatively, using linear response theory 
In this model, two localized quantum states correspond- 
ing to the redox sites are taken. This two-level system 
(TLS) is fully specified by a tunnel splitting HA (which 
is twice the electronic coupling between the two redox 
sites), and by a bias he (which corresponds to the energy 
difference between the sites). Of crucial importance is 
then the coupling to the "environmental" modes (bath), 
which is modeled as an infinite set of harmonic oscilla- 
tors with a suitable (continuous) spectral density J{lo) 
describing their frequency-resolved coupling strength to 
the electron. For a given ET system, one can use classical 
Molecular Dynamics simulations to find the appropriate 
J(lu) [El, but below we shall consider a specific class of 
spectral densities that describes many ET processes quite 
accurately. Two important parameters following from a 
given spectral density are the reorganization energy HA 
(which describes the overall coupling strength) and the 
typical bath frequency uj c corresponding to a maximum 
in J(lo). This model has been studied in depth in the 
field of dissipative quantum mechanics f7|-|lT|| , mostly in 
the so-called scaling limit for Ohmic spectral densities 
(see below), but its application to ET reactions is mainly 
located in a different parameter regime. 

Certain limiting cases of ET theory are rather well un- 
derstood. The primary object of interest is the thermal 
rate constant kth, which is useful if ET dynamics is cor- 
rectly described by exponential decay. First, for very 
small A, one can use perturbation theory in A to ob- 
tain the celebrated golden rule rate JI^15|. This is the 
nonadiabatic limit. In the opposite limit of large A, the 
adiabatic limit, it is more useful to think in terms of adia- 
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batic free energy surfaces formed by linear superpositions 
of the two localized free energies. In that case, the bath 
essentially is very slow ( "classical" ) compared to the elec- 
tronic motion, and one can develop a simple theoretical 
description again jl6| Q. Notably, the ET dynamics in 
this regime is often not simply an exponential relaxation 
but exhibits oscillatory or multi-exponential behavior at 
low temperatures |H^0|. At sufficiently high tempera- 
tures, the famous classical Marcus theory or semi- 
classical generalizations thereof S, apply. For a path- 
integral approach to Marcus theory, see also Ref. 
Marcus theory also makes predictions about the rate con- 
stant in the crossover regime between nonadiabatic and 
adiabatic ET. It should be stressed that the spin-boson 
model contains all these results as limiting cases in the 
full parameter space. 

Notably, the spin-boson problem cannot be solved an- 
alytically without imposing approximations. One impor- 
tant exception to this rule is the case of an Ohmic spectral 
density, see Eq. ([l7|) below, with the restrictions: 

A/w c «l, k B T/fiuj c < 1 <^ scaling limit . (1) 

In this scaling limit, sophisticated conformal field theory 
techniques are able to yield the complete non-equilibrium 
solution for the dynamics under arbitrary initial prepara- 
tion p^| . Besides this solution, whose parameter regime 
is of little interest to ET theory, if one is interested in 
analytical progress, one has several options for approxi- 
mations. One possibility is to work perturbatively in the 
tunnel splitting A. This leads to the noninteracting blip 
approximation (NIBA) [p3|-p5|, which is basically equiv- 
alent to the golden rule formula for the rate constant. 
Other approaches are based on variational treatments 
p6| or work perturbatively in the system-bath coupling 
(Redfield approach) . As typical reorganization ener- 
gies are quite sizeable, a perhaps more promising route is 
to use the powerful computer simulation methods avail- 
able by now. 

Let us then briefly summarize the numerical meth- 
ods that have been applied to a study of the spin- 
boson dynamics. An approximate approach which 
has become quite popular recently employs a mixed 
quantum/classical (or quantum/semiclassical) descrip- 
tion |2§|-|37][. However, such studies are typically re- 
stricted to a few bath modes. Unless one has a sit- 
uation with only a few strongly coupled bath modes, 
the more typical situation of condensed-phase ET re- 
actions is difficult to describe using such methods. In 
addition, it is hard to estimate the errors made un- 
der these approximations. Similar arguments apply to 
time-dependent Hartree methods (basis set approach) 
p8| and to memory-equation approaches J3£|[40]]. The 
latter method, as well as the recently proposed stochas- 
tic Langevin-type simulation flllp^ ], dynamical Wilson 
renormalization group |43||44| ] and the real-time renor- 
malization group Ha], basically appear to be restricted 



to the Markovian limit realized for an Ohmic bath with 
lu c — > oo, yet most situations of interest for ET are lo- 
cated outside the Markovian regime. The methods dis- 
cussed so far are therefore very powerful in certain pa- 
rameter regions of the spin-boson model (in fact much 
more powerful than our method employed below), but 
may fail badly in others. 

A method free from such restrictions is provided by 
Monte Carlo simulations p6[ , in particular the powerful 
path-integral Monte Carlo (PIMC) approach. In princi- 
ple, this approach is numerically exact throughout the 
full parameter space, i.e. for arbitrary spectral density, 
tunnel coupling, temperature and/or bias. However, as 
one is really interested in the ET dynamics in real time, a 
direct application of quantum Monte Carlo (QMC) meth- 
ods is plagued by the infamous sign problem [147). One 
alternative is to carry out an imaginary-time simulation, 
which is intrinsically free of the sign problem, and then 
analytically continue the data to real time. This program 
has been carried out for the spin-boson problem, cither 
using Pade approximants Jl8],fl9| or Maximum Entropy 
methods J50| ] in order to perform the analytical continu- 
ation. The latter procedure is mathematically ill-defined 
and quite troublesome to carry out. In addition, by us- 
ing this procedure, one is inevitably restricted to equi- 
librium situations. However, in ET reactions one often 
has a true non-equilibrium preparation, e.g. one photo- 
excites an electron to the donor state at time t = and 
then follows the time evolution of the occupation prob- 
ability. Such a question has to be answered by a true 
non-equilibrium real-time calculation. 

Such a calculation is presented in our study. The price 
to pay for working in real time and still keeping exactness 
and applicability throughout the full parameter space is 
however quite high. The main challenge is due to the sign 
problem, which implies that one is restricted in the real- 
time range up to which stable calculations can be done. 
The sign problem arises from the destructive interfer- 
ence of a large number of quantum paths at long times, 
leading to a very small signal-to-noise ratio and thus ren- 
dering the simulations numerically unstable. Early real- 
time PIMC attempts for computing the spin-boson dy- 
namics J5l]-|53| were based on the stationary-phase Monte 
Carlo approach |54|-|5(| or variations thereof. Such filter- 
ing methods try to keep the MC trajectory close to the 
important stationary-phase regions of path space. While 
this is able to somewhat relieve the sign problem, a more 
powerful method is available by now. This "multilevel 
blocking" (MLB) algorithm [|57|-|60|| has been proposed 
and successfully used to relieve both the fermionic and 
the dynamical sign problem. MLB represents a system- 
atic recursive implementation of a simple blocking strat- 
egy. Roughly speaking, the blocking strategy is equiva- 
lent to a naive filtering, and one can think of MLB as a 
systematic and optimized way of implementing the filter- 
ing idea. Here we shall employ a recently proposed MLB 
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variant which can account for effective actions BiJ of 
the type encountered in PIMC for the spin-boson model. 
This approach is superior to previous real-time PIMC 
methods, and allows to compute the spin-boson dynam- 
ics up to timescales of practical interest without approxi- 
mations and for the full parameter regime relevant to ET 
reactions. 

The remainder of this paper is as follows. In Sec. we 
briefly describe the spin-boson model and the dynami- 
cal quantities of interest. Next we describe our QMC 



method, see Sec. III. Since the method has been exposed 
in great detail in Ref. |5S[], focusing on exactly the same 
model, we largely refrain from repeating ourselves here. 
Technical details of interest to experts are given in the 
Appendix. In Section IV we then present our numerical 
results for the dynamical quantities and the rate kth for 
various regimes, including the important crossover region 
between nonadiabatic and adiabatic ET. The most dif- 
ficult regime is at low temperatures, a regime that can 
now be treated by PIMC. In presenting the numerical re- 
sults, we restrict ourselves to the symmetric (unbiased) 
case e = 0, leaving the biased system for future studies. 
Finally, some conclusions are drawn in Sec. 



II. ELECTRON TRANSFER THEORY 

The spin-boson model (dissipative two-state system) is 
defined by the Hamiltonian 

H = H + H I + H B = -^<?x + 
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where the free TLS is described by Hq, containing the 
tunnel splitting hA and the energy gap he between the 
two localized electronic states. Here a x and a z denote 
the standard Pauli matrices acting in TLS space, where 
the |+) (|— )) state refers to the donor (acceptor). The 
environmental modes are modeled by an infinite collec- 
tion of harmonic oscillators, Hb, which bilinearly couple 
to the position of the electron (Hi). For the great ma- 
jority of ET systems, this model provides a reasonably 
accurate description of reality, as is discussed and moti- 
vated at length in Rcfs. ]5| |7|,|l5|. Within this model, the 
bath parameters affect ET only via the spectral density 
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which effectively becomes a continuous function of to for 
a condensed-phase environment and determines all bath 
correlation functions that are relevant for the ET dynam- 
ics. For instance, the complex-valued bath autocorrela- 
tion function for complex time z = t—ir is for (3 = 1/ksT 
given by 



L(z)= f^M^M^M- (4) 
w J 7T y ' sinh(^/3w/2) v ' 

Similarly, the important reorganization energy hA of 
Marcus theory 0,^] is defined by 



A 



dbJ 
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which is the only important bath quantity in the classi- 
cal (high-temperature) limit, as is apparent from classical 
Marcus theory. Details of J(u>) matter only at low tem- 
peratures. 

In the following we focus on the two dynamical prop- 
erties of primary relevance to ET dynamics. The time- 
dependent occupation probability is defined as 



P(t) = (a z (t)) = (e 



_ iiHt/h -iHt/h 



(6) 



and gives the difference in the occupation probabilities 
of the donor and the acceptor state with the electron 
initially held fixed on the donor. This quantity then 
directly probes ET dynamics after the non-equilibrium 
initial preparation corresponding to the initial density 
matrix 



W(0) = \+}(+\e- fi{HB+ ^ 



(7) 



Here fi is the dipolc moment of the electron, and £ de- 
notes the dynamical polarization of the bath Q, with 
the electron in the donor state and the bath equilibrated 
with respect to it. By comparing with Eq. (|J), we see 
that \xE = — J2 a C a X a . As pointed out in Ref. Q, this 
"standard preparation" often used in ET experiments 
is unfavorable for coherent (oscillatory) dynamics when 
compared to other initial preparations. Therefore it is es- 
pecially suitable for a thermal transfer rate description. 

In addition, we calculate the complex-valued correla- 
tion function 



C(t) = (a z (0)a z (t)) fi 
= Z-Hr{ 



e-P H o z e iHt l h o z e- lHt l h 



(8) 



with Z = tr exp(—[3H). This quantity probes the dy- 
namics under an equilibrium preparation and, for e = 0, 
differs from P(t) only through the initial preparation. 
At high temperatures, one generally expects that initial 
preparation effects show little influence on the dynamics, 
and therefore a unique rate constant kth should exist. 
The inverse rate would then be the timescale governing 
the relaxation of both P(t) and C(t). However, in the 
quantum regime one needs to be more careful, and the 
specific experiment under study will determine which is 
the relevant dynamical quantity. 

Assuming for the moment that there is a unique rate 
constant k t h, P{t) should follow a simple exponential de- 
cay for times t large compared to some transient molecu- 
lar timescale T tran s of the order l/w c , yet small compared 
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to the relaxation timescale r ro i ax » l/kth- Calculating 
P(t) for T trans «C t <C r rc i ax would then allow to extract 
the temperature-dependent thermal transfer rate. The 
thermal transfer rate (if it exists) can also be obtained 
via a time-dependent function defined under an equilib- 
rium preparation (Hi , 



*/(*) 



-Imtr [e-P B h A (0)h A (t)] 



(9) 



where h A = |+}(+| = (l + 0z)/2 denotes a projector onto 
the donor state and Z A = tr{e~@ H h A } its partition sum. 
The thermal rate constant then follows from Eq. (||) as 
the plateau value for times t where fe/(t) forms a plateau. 
According to our above discussion, this should happen for 
Ttrans *C tp] <C T ro i a x where £ = ipi is a suitable plateau 
time. The index / in Eq. (g) denotes a forward rate, 
describing the directed transfer rate from the donor to 
the acceptor site. The total rate kth probed by P(t) is 
the sum of the forward and the backward rate kf and 
fcft, which in turn are connected via a standard detailed- 
balance relation, 



k b (T, e) = kf(T, e) exp(-?ie/fc B T) . 



(10) 



Notably, if T trarls and r re i ax are not well separated, no 
transfer rate can be defined, and a simple rate formal- 
ism breaks down. Such a breakdown of the rate descrip- 
tion could happen, for instance, for very fast reactions. 
Equation (|9|) thus not only yields a convenient way to 
compute the thermal transfer rate, but also a means to 
decide upon the validity of the rate picture of ET at all. 
It is important to stress that kt(t) does not refer to a 
time-dependent rate reflecting the dynamics of P(t) but 
is only an auxiliary function allowing to determine the 
time-independent rate kth when it exists. Equation (||) 
can be expressed in terms of C(t), 
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1 Im C(t) 



(11) 



with (a z )/3 = Z-Hx{e- 0H a z }. The fact that P(t) will 
follow an exponential decay only after some time T trans 
is made explicit by noting that dP/dt(t = 0) = due to 
da z (t) /dt ~ a y (t) . Our main goal is therefore to calculate 
the ET dynamics up to real times that are significantly 
larger than r trans . 

Before turning to the numerical computations, let us 
briefly summarize some simple limiting cases, where an- 
alytical results are available. In the nonadiabatic limit, 
Eq. (^) can be evaluated for arbitrary temperatures using 
Fermi's golden rule, 



k GR 
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dt exp[iet — Q(t)] 



(12) 



with the twice-integrated bath autocorrelation function: 



d 2 Q(z) /dz 2 = L(z) with Q(0) = . 



(13) 



For special choices of J(ui), it is possible to evaluate the 
remaining integral explicitly ]l5| . In the adiabatic limit, 
the bath is very slow and thus behaves classically. Es- 
sentially, it then represents a static random field obeying 
Gaussian statistics, which acts as an additional fluctu- 
ating energy gap between the electronic states. The ET 
dynamics then follows from the free TLS dynamics by 
including a fluctuating bias || . The inevitable failure of 
this approach at long times can in principle be relieved by 
applying appropriate Bloch-type equations Jig]. Third, 
for high temperatures, classical Marcus theory Eq] can be 
recovered from the spin-boson description plj and pre- 
dicts the forward rate 



k c l = 



A 2 



/ 4 + 7rA 2 /(Aw r ) V Ak B T 



-F*(e)/k B T 



(14) 



with the activation free energy barrier given by the cele- 
brated Marcus parabola, 



F*(e) = h(e- A) 2 /4A, 



(15) 



and a solvent frequency scale uj r . For the Ohmic spectral 
density (|17]), this frequency was computed approximately 
in Ref. ElTL with the result 
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(16) 



Marcus theory yields a classical rate constant covering 
the full crossover from nonadiabatic to adiabatic ET. 
Equation (|lj) coincides with the high-temperature limits 
of the nonadiabatic rate (|l|) and of the adiabatic formu- 
lation for A -C lo c and A 3> u> c , respectively. In the adia- 
batic limit, the rate is independent of A. ^From Eq. JT^ ), 
the crossover regime is expected for A/uj c \J A/u> c , at 
least for high temperatures. 

In all simulations reported in this paper, we consider 
an Ohmic spectral density with exponential cutoff M , 



(17) 



with dimcnsionless damping strength a and a cutoff fre- 
quency u> c . For many polar solvents, it is appropriate to 
choose uo c according to some intermediate bath frequency 
]l5[ . Instead of a, it is often more convenient to measure 
the coupling strength to the bath in terms of the reorga- 
nization energy (^). For this spectral density, one finds 
A = 2auj c . In addition, the correlation function ( |l3| ) can 
be given in closed form Q , 



Q(z) = 2a 



-In 



ln(l + ioj c z) 
r(fi + iz/h(3)T(n - iz/h/3) 



(18) 
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with the gamma function T(z) and the abbreviation 

n = l + k B T/hu c . (19) 

Note that the second contribution in Q(z) in Eq. JTs| ) 
vanishes for T = 0. 

Let us also summarize at this point some of the ana- 
lytically known results for this class of spectral densities. 
First, the complete temperature dependence of the nona- 
diabatic golden rule rate ( |l2| ) is known for a = 1/2 and 
a = 1. For e = 0, the result is f0~5|] 



kf R (a = 1/2) 
kf R (a = 1) = 



ttA 2 £(20 - 1) 

Wc 2 2n r 2 («) ' 
ttA 2 r 4 (2Q - 1) 

2n/3w 2 T 4 {n)T{Afl - 2) 



(20) 
(21) 



Second, the complete dynamics for a = 1/2 in the scal- 
ing limit defined by Eq. (|l|) is known. For instance, the 
quantity P(t) obeys a simple exponential relaxation for 
e = and Lo c t 3> 1, 



(22) 



where the relaxation rate is twice the low-temperature 
limit of Eq. @. Notably, for a ^ 1/2 but close to 1/2, 
the dynamics is much more complicated, with algebraic 
decay and oscillatory behaviors Q. These features are 
due to electronic coherence effects, which one generally 
does not expect to survive in ET reactions. This is also 
reflected in the behavior of C(t) for a = 1/2, which ex- 
hibits long-time algebraic tails. Therefore one cannot ex- 
pect that the time-dependent function kf(t) will exhibit 
plateau behavior, although it is possible to define a rate 
for the decay of P(t). This reflects the importance of the 
initial preparation in that case. 



III. SIMULATION METHOD 

Next we turn to a computational method able to yield 
the rate constant. A great simplification arises due to the 
harmonic nature of the bath and the bilinear coupling in 
Eq. (|J), which allows to perform the trace over the bath 
degrees of freedom in Eqs. (|^) and (|J) in an exact manner. 
Adopting a path-integral approach, one obtains, e.g. for 

c{t) 0, 



C{t) = z~ 



Va a{0)a{t) exp <{ r^oM - $M } . (23) 



Here the path integration runs over paths o[z) for the 
discrete "spin" variable a = ±1 corresponding to er 2 , 
with the complex time z following the Kadanoff-Baym 
contour 7 depicted in Fig. |l|. Furthermore, So[a] denotes 
the total action of the free TLS, and the influence of the 



traced-out bath is encoded in the influence functional 
<f>[cr]. In terms of the bath autocorrelation function ([|), 
it reads 



$[cr] = - I dz dz'a{z)L{z - z')a(z') , 

4 Jz'<z 



(24) 



where time integrations are ordered along the contour 7. 
The influence functional introduces long-ranged interac- 
tions along the real and imaginary time axes, rendering 
the evaluation of the remaining path integral in Eq. ( |23| ) 
a difficult task. 

Quantum-mechanical expectation values like Eq. ( [23] ) 
can be calculated in a numerically exact way by PIMC 
simulations. To make Eq. (|23| ) accessible to PIMC, the 
path integral is taken in its discretized form with N 
discretization steps. Following the scheme described in 
Refs. |53|,^{|, q uniformly spaced points are used for each 
of the real-time paths z : — > t and z : t —> 0, and r 
points for the imaginary-time path z : — ► —ihf3, with 
discretization steps 



t/q, !<3<q 
-t/q, q+l<j<2q (25) 
-ih/3/r, 2q+l <j <2q + r , 



leading to N = 2q + r discrete time points Zi = J2]=i A? 1 
Using Oi — a(zi) = ±1, Eq. (|23j) can be written as 



C(t) = Z 



_1 E 



N 



i=l 



-*[{^}] 



(26) 



with Z = } P[{ a j}} an d cjv+i = 01 due to the cyclic 
nature of the trace. The sum runs over all realizations of 
the discretized spin path {<Jj} — {ci, . . . , and 



K(ai,a i+ i) = (a l+1 \exp(~iAiH /h)\ai) 



(27) 



denotes the short-time propagator of the free TLS, which 
is a known 2x2 matrix. Note that due to {a(0)a(t))p = 
{a(t')a(t' + t))p, one single MC trajectory can be used 
to compute C(ti) for all times t\ < t. Furthermore, 
exploiting (a(0)a(tk)) p — (a(tk)a(0))% in a similar way 
improves the MC statistics. A discretized version of the 
influence functional is given by 



*[fo}] 



1 ^ 



j,k=l 



(28) 



with the complex-valued influence matrix |p3| 

L jk = L kj = Q(zj -z k + (Aj + A fc _i)/2) 
+Q(zj -z k + (-Aj_i - A fc )/2) 
-Q(zj -z k + (-Aj_i + A fc+ i)/2) 
-Q(zj -z k + (Aj - A fe )/2) for j>k, 

L Jj - = 2Q((A i _ 1 + A J -)/2), (29) 
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where Q{z) is given in Eq. To make further 

progress, and also to correct a misprint in Ref. ]55|], we 
now denote the spins a2q+2-j for 1 < j < q + 1 (residing 
on the backward real-time branch of 7) by a'j , and, sim- 
ilarly, the imaginary-time spins (J2q+m for 2 < m < r by 
a m . In the next step, employ the coordinate transforma- 
tion 

Vj = + v'j) , = ~ v'j) ■ 

These are "classical" and "quantum" variables for elec- 
tronic motion j53j. Exploiting the symmetry relations 
Q{z - ih(3) = Q(-z) and Q{z*) — Q*(—z), one finally 
arrives at 



r ^ q 

rn > n — 2 



j>k=l j = l m=2 

a' r 

m,2g+l) , 

m=2 

with the matrices 



(30) 



Yrnn — ^p2q+m,2q+n j 



Xjfe = Im L 



Ajk — Re Ljk 
1 



where 1 < j, k < q and 2 < m,n < r. Equation ( J3C| ) 
together with the contribution from the TLS propaga- 
tor ( p7| ) constitutes a discretized form of the total action 
which is useful for PIMC simulations. Finally, to obtain 
P(t) instead of C(t), the imaginary-time spin contribu- 
tions in Eq. ([30j) and in the free propagator have to be 
neglected, with a\ = a\ = +1 reflecting the initial non- 
equilibrium preparation. Instead of sampling an entirely 
new MC trajectory according to this new weight, P(t) 
can also be calculated from the same MC trajectory as 
C(t) by including a correction factor which accommo- 
dates the changes in the corresponding MC weight. We 
refer the interested reader to the Appendix for compu- 
tational details. Although the latter approach exhibits 
poorer statistics, it still produces satisfying results in 
most cases. We therefore use it for the simultaneous cal- 
culation of kf{t) and P{t), without significant increase in 
computing time. 

Unfortunately, this method suffers from the notorious 
sign problem. It stems from interference due the com- 
plex phase factors assigned to different spin paths {(Jj}, 
resulting in a small signal-to-noise ratio of the stochastic 
averaging procedure. The exponential increase of inter- 
fering paths with system size N is reflected in an ex- 
ponential decrease of the signal-to-noise ratio with the 
maximum real time under study. Here we employ the 



multilevel (MLB) approach in a version suitable to deal 
with long-ranged interactions along the (complex) time 
contour, which is capable of relieving the sign problem 
without introducing approximations. The algorithm is 
described in detail in Ref. |5{| , and we only give the basic 
idea in the following paragraph. For the expert reader, 
some improvements over Ref. [ |59| are summarized in the 
Appendix. 

In short, the MLB algorithm is based on the simple 
observation that the sign problem can be considerably 
reduced by sampling blocks instead of single paths. A 
block is here built from a small set of paths. Provided 
the blocks are chosen small enough, stochastic averages 
within a block do not suffer from the sign problem. Sam- 
pling blocks instead of single paths can in fact be shown 
to always result in a better signal-to-noise ratio |5(| . Re- 
peating this scheme iteratively by establishing a hierar- 
chy ("levels") of blocks, the exponential severity of the 
sign problem with respect to the maximum real time t 
can be changed into an only algebraic one j57|-|6Ct| . 

Single spin flips and, on the real-time axis, double spin 
flips, with simultaneous flips of the forward and back- 
ward spins <Ji and a[, were used to propagate the MC 
trajectory. Furthermore, since the creation of short blips 
(connected real-time intervals where & = 0) is often en- 
ergetically more expensive than the creation of more ex- 
tended ones, we also included flips of complete blocks 
of subsequent spins to ensure ergodicity. When working 
with MLB, these blocks must not extend over the border- 
line separating different levels. However, since the results 
do not deviate from those obtained with local spin flips 
alone, we conclude that ergodicity poses no problem for 
our approach. 



IV. NUMERICAL RESULTS 

In this section, we present numerical results for the 
thermal transfer rate kth and the occupation probability 
P(t) of an unbiased system, e = 0. Of course, then kth 
is just twice the forward rate kf. Accordingly, kth wa s 
extracted from the plateau value of 



k(t) - 2k f (t) . 



(31) 



Furthermore, the right-hand side of Eq. ( |11|) reduces to 
lm{C(t)} /h(3 since (a z )p vanishes for e = 0. Our code 
was carefully tested for (a) a = (free TLS), (b) a = 1/2 
within the scaling limit (^), and (c) for a = 5, uj c — 
0.5 A, T = 2hA/kB- The numerical results accurately 
reproduce the available analytical solutions [for (a) and 
(b)] and the numerical results of Ref. j3^] in case (c). For 
more tests, see Sec. [VA. 



We are then interested in mapping out the ET dynam- 
ics phase diagram, which is a function of three parame- 
ters. These are (1) the temperature T, (2) the reorga- 
nization energy A, where for the bath spectral density 
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( |l7j ) one has A = 2au c with a damping parameter a 
and a cutoff (dominant) frequency u> c , and (3) the adia- 
baticity parameter A/cu c with the tunnel splitting A. We 
shall systematically tune these parameters in order to ex- 
plore the crossover from nonadiabatic to adiabatic ET. 
To keep the relevant parameter space manageable, we 
restrict ourselves to the reorganization energy A = 10 A, 
but also give some results for A = 50 A |6l|] . Unless noted 
otherwise, below A = 10A has been taken. We are then 
left with only two free parameters, namely temperature 
and the adiabaticity A/lo c . By systematically increasing 
the latter, one can go from the nonadiabatic into the adi- 
abatic region of the phase diagram. Of course, we use the 
term "phase diagram" not in the sense of distinct regions 
separated by phase transitions. Nevertheless, there are 
regimes where ET dynamics proceeds qualitatively dif- 
ferent, which are separated by a rather sharp crossover. 

The thermal transfer rate was numerically obtained in 
two different ways. If k(t) showed a plateau for some 
plateau time t p \ > r trans , the thermal transfer rate kth 
was taken as the plateau value k(t p \). A second and dis- 
tinct possibility is to fit P(t) to an exponential decay 
~ exp(— k t ht) after the initial transient, t > Tt ran s- In 
fact, we will encounter examples where a rate from k(t) 
does apparently not exist, or at least requires a modi- 
fied extraction procedure, while P(t) still displays a well- 
defined exponential decay on the timescale under study. 

A. Special cases 

To make contact to analytical theory, let us first de- 
scribe results for the special damping strengths a = 1/2 
and a = 1. This provides also another check for the nu- 
merical MLB-PIMC scheme. In addition, from our data 
for a = 1/2, we can give fairly accurate estimates for the 
validity of the scaling regime, namely 



w c /A>6 hoj c /k B T>50 



(32) 



Let us start with a = 1/2 and lo c — 50A, corresponding 
to the large reorganization energy A = 50 A. Here an 
almost strict exponential decay of P(t) was found at all 
temperatures, with the rate nicely matching the analyti- 
cal prediction (20), see Fig. ^. This agreement is of course 
not surprising since lo c is much larger than A. For P(t), 
transient short-time dynamics reflecting dP/dt(0) = 
takes place on a significantly shorter timescale than for 
k{t), where t p \ increases for lower temperatures. For 
T <TiA/kB, k(t) fails to exhibit a plateau while P(t) 
still allows for a useful rate description. This shows the 
increasing significance of electronic coherence in this pa- 
rameter regime . 

This point becomes even more obvious for A = 10A, 
where A/u c = 0.1. Now k(t) fails to exhibit a plateau at 
all temperatures studied, 0.2 < ksT/fiA < 10. How- 
ever, the rate obtained from P(t) still closely follows 



the nonadiabatic prediction (pp|), see Fig. | In addi- 
tion, for a = 1/2 and T = hA/ke, the dependence on 
the adiabaticity parameter was studied within the range 
10 < uj c / A < 100. Again the nonadiabatic prediction is 
nicely reproduced, proving its usefulness in this regime, 
see Fig. 0. The excellent agreement of Eq. (^TJ) with our 
results allows for the quantitative identification of the 
scaling regime. The rate ( p0| ) reaches the scaling value 
7rA 2 /2w c within a relative error of r %, if 



huj c 
k^f 



423 
r 



(33) 



Our data indicate that this criterion also holds for the 
numerically exact ET rate for (at least) 0.1 < r < 10, 
0.01 < A/u c < 0.2, and 0.05 < k B T/T%A < 20. Note 
that this also agrees with Eq. ([32]) for r ~ 8.5 %. 

The second special case, to which we now turn, is de- 
fined by a — 1. Again starting with A = 50 A, the ob- 
tained rate closely follows the nonadiabatic prediction 
( pl| ) for all temperatures investigated, 0.2 < fc^T '/TiA < 
20, see Fig. whereas the classical (high-temperature) 
Marcus rate (|l4j) fails to give an accurate picture ex- 
cept for T > AQfiA/ks = 1.6Tiuj c /kB- The absence of 
electronic coherence and the nice separation of transient 
and relaxation timescales is reflected in the fact that k(t) 
reaches a plateau within the whole temperature range. 
The corresponding thermal rate constant is identical to 
the rate extracted from the exponential decay of P(t > 
Ttrans)- Moreover, in contrast to a = 1/2, this exponen- 
tial decay is preceded by a steep decrease of the popula- 
tion for T < 5hA/ks on timescales t< T\, TWS ~ 0.5A -1 . 
The timescale rt rans (approximately) coincides for both 
P(t) and k(t), see Fig. |, and remains fairly constant with 
temperature. The a — 1 rate for given other parameters 
is smaller than for a = 1/2. Since with our conventions, 
a = 1/2 corresponds to a "more nonadiabatic" regime, 
this finding is not surprising. Recrossing events, which 
are irrelevant in the nonadiabatic limit, will generally de- 
crease the rate. 

For A = 10A, we find a similar scenario, with the Mar- 
cus rate failing except for T > lOftA/fc^ = 2fiuj c /kB- 
However, k(t) does not exhibit a plateau, since ET dy- 
namics is now very fast and hence Tt rans <C T re i ax does 
not hold anymore. Remarkably, P{t > T trans ) still shows 
exponential decay with the nonadiabatic rate (pl|), see 
Fig. |j. Note that here the adiabaticity parameter is al- 
ready as large as A/lj c = 0.2. 



B. Crossover regime 

Entering the crossover regime between nonadiabatic 
and adiabatic ET (taking A = 10A), let us start with 
A/oj c = 1. In this case, the extraction of the ther- 
mal transfer rate from kit) fails for high temperatures, 
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T > j 2hA/kB- After an initial peak, k(t) decreases al- 
most linearly instead of exhibiting a plateau. The high- 
temperature absence of a plateau may be rationalized 
by noting that the classical activation barrier ( plj ) corre- 
sponds to T = 2.5hA/ks- For higher temperatures, ET 
is practically activationless, so that transient and relax- 
ation timescale are not really separated anymore. How- 
ever, after a rapid initial transient decay, P(t > Tt ran s) 
decays exponentially with the rate kth- Interestingly, this 
rate coincides with the value of k{t) at the transition from 
the initial peak to the linear decrease, see Fig. ^, and gov- 
erns the behavior of P(t) for rt rans < t < T tr ans + O-Sfc^, 
(note the increase of the relevant time interval for low 
temperatures.) The temperature dependence of k t h is 
well described by the Marcus formula (jli|), see Fig. 0. 
While the Marcus rate together with Eq. ( jltf ) yields a 
quite accurate description, the nonadiabatic rate Jl2] ) 
overestimates the true rate by a factor rj 1.2. How- 
ever, overall features (like the temperature at which the 
maximum of the rate occurs) are still reproduced. This 
overestimate of the rate can be rationalized in terms of 
recrossing events between acceptor and donor, which are 
higher-order contributions in A. 

For low temperatures, T < 0.2hA/kB, we have encoun- 
tered a rather interesting and apparently quite general 
phenomenon regarding the ET dynamics in the crossover 
regime. In fact, in this low-temperature crossover regime 
between nonadiabatic and adiabatic ET, the thermody- 
namic rate becomes so small that almost no change in 
P(t > T trans ) is observed anymore. At the same time, the 
initial decay of P{t) is ended by a local minimum deep- 
ening with decreasing temperature, see Fig. |8|. Therefore 
the relevant ET dynamics solely happens during the ini- 
tial (transient) stage, which lasts for a few A -1 , while the 
subsequent approach to equilibrium takes place on an ex- 
tremely slow timescale. This freezing of the ET dynam- 
ics was never observed in the nonadiabatic regime, but 
prevails both in the crossover and the adiabatic regime. 
Below we refer to this novel behavior as "transient ET," 
in contrast to the conventional limiting cases of adiabatic 
and nonadiabatic ET. 

Next we increase the value of A/oj c to 2. The nu- 
merical results in Fig. |^ reveal a similar picture. Again 
k(t) fails to exhibit a clear plateau for T ^hA/ks, but 
kth can be extracted from P(t). The Marcus formula 
( fli] ) captures the temperature dependence of kth quite 
well, see Fig. ^, while the nonadiabatic prediction still 
reproduces qualitative features, now overestimating the 
rate by a factor « 1.4. Adiabatic dynamics jl6| already 
describes the transient dynamics, underlining adiabatic 
tendencies in the bath, but no oscillations are observed 
(within error bars). For T <0.6hA/kB, P(t) develops 
a local minimum-maximum pair in the short-time do- 
main t< A" 1 , whose amplitude increases as T — > 0. 
For T < 0.3hA/ks, we again observe transient ET as de- 
scribed above. Apparently, ET dynamics becomes quite 



complex in this regime. 

C. Adiabatic limit 

The transition into the adiabatic regime was stud- 
ied by increasing A/u> c up to 20 for two temperatures, 
T = 10hA/k B and hA/k B - Both P(t) and k(t) are 
expected to eventually follow the adiabatic prediction 
]l6| for large A/u> c , showing damped oscillations with 
decreasing amplitude for lower temperatures. These os- 
cillations are not due to electronic coherence but a con- 
sequence of "nuclear" (vibrational) coherence. Although 
adiabatic theory correctly gives dP/dt(0) — 0, it is well- 
known not to reproduce the correct long-time dissipa- 
tive behavior 0. Nevertheless, for t <C r rc i ax , conven- 
tional wisdom expects that it should still provide a useful 
approximation. The short-time dynamics of both P(t) 
and k(t) indeed closely follows adiabatic theory. For 
T = HA/ks, a plateau in k(t) is found, yielding a trans- 
fer rate kth that decreases with increasing A/u> c , until 
oscillations eventually occur at A/oj c >5. Since (within 
error bars) P(t) shows no oscillations even for A/tv c = 5, 
we can assign a rate to ET under a nonequilibrium initial 
preparation. Only for larger A/u> c , evidence for oscilla- 
tions in P(t) can be observed, see also Ref. p9| . 

Several important differences to adiabatic theory ap- 
pear in our data, see Figs. |l(] and O. With regard 
to k(t), the observed oscillations have a smaller ampli- 
tude, their mean value is larger, and the oscillation fre- 
quency is higher. Adiabatic theory assumes a static bath, 
which under a dynamical description implies incorrect 
long-time behavior and overly pronounced oscillations. 
This qualitatively explains the first two differences. The 
increase of the oscillation frequency is then due to the 
larger energy difference between the two adiabatic po- 
tential surfaces away from the barrier region, since the 
region away from the barrier top is probed extensively 
by a dynamical bath. Turning next to P(t), adiabatic 
theory fails to reproduce the observed steep descend at 
the end of the transient dynamics. This can be explained 
in terms of the bath distribution. In adiabatic theory, the 
distribution remains centered around the donor state, see 
Eq. (0) , so that there is a low probability of reaching the 
transition region, and hence too slow ET dynamics is pre- 
dicted. The improved adiabatic description |]l6| captures 
the damping of the oscillation amplitude and the increase 
of the frequency much more accurately. However, since 
it is still based on an average over a static bath distri- 
bution, it again misses the steep descend of P(t) and 
the correct mean value for k(t) at times t > ttrans- For 
increasing A/u) c , adiabatic theory then provides a sys- 
tematically better description, although substantial dif- 
ferences remain even at A/u> c — 20. For T — lOhA/ke, 
we get a similar picture, see Fig. [w]. However, as ex- 
pected from our above discussion, k(t) does not exhibit a 
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plateau. The transition from incoherent to coherent (os- 
cillatory) behavior in P(t) takes place between A/oj c = 5 
and 10. Remarkably, true adiabatic dynamics is realized 
only for about half an oscillation period, this being virtu- 
ally independent of A/ui c . This shows that the long-time 
relaxation process not captured by adiabatic theory sets 
in much earlier than thought previously. 



D. Validity of rate description 

Having identified the threshold to the adiabatic regime 
(at least for A = 10A and k B T/hA = 1,10), we next 
comment on the validity of a rate description for ET. A 
rate description for P(t) turns out to be appropriate for 
A/uj c <5, where accurate fits of the numerically exact 
rates to analytical theory are possible. For small A/oj c , 
the rate is given by the nonadiabatic prediction (12), 
while in the crossover regime, the Marcus rate ( |l4| ) is su- 
perior. The Marcus formula assumes a classical bath, and 
hence becomes more appropriate for large A/u> c , where 
the bath is very slow and hence classical. Remarkably, 
the Marcus formula can be made to work even outside 
the true classical (high-temperature) regime, see Fig. |lj. 
In the high-temperature limit, the solvent frequency tu r 
appearing in Eq. ([u]) is given by Eq. (|g|). For lower 
(but not too low) temperatures, the Marcus formula still 
provides a good estimate for the exact rate within the 
regime stated above, provided uj r is taken in the form 



Y 



(x) 



1{T) 



(34) 



In the high-temperature limit, our data indicate that 
Eq. ( |34"| ) reduces to Eq. (Hlj), implying q(T) — > for 
high temperatures. At low temperatures, however, a pos- 

Of course, at 
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itive value for q(T) is found, see Fig. 
sufficiently low temperatures, the quantum-mechanical 
"transient regime" mentioned above is entered, where 
the Marcus formula breaks down. Nevertheless, Eq. (|34| ) 
represents a noticeable improvement over Eq. ( |l6| ) for a 
significant range of intermediate temperatures. We term 
this low-temperature region, where classical Marcus the- 
ory along with Eq. ( |34]) becomes accurate, "extended 
Marcus regime." That such a regime could indeed ex- 
ist may be rationalized by noting that for A/u> c >l, the 
bath behaves already rather classically, yet quantum fluc- 
tuations significantly renormalize the solvent frequency 
scale u> r at low-to-intermediate temperatures. It would 
clearly be interesting to provide an analytical derivation 
for Eq. @. 

Notably, for A/w c >, 10, no rate description is possible 
anymore, neither based on P(t) nor on k(t), since ET dy- 
namics is dominated by vibrational coherence. The rate 
formalism also becomes useless at very low temperatures, 
T <0.2HA/kB, within the crossover regime. As elabo- 
rated above, in this "transient regime," kth is almost zero 



such that long-time relaxational ET is essentially frozen 
in. The important ET dynamics then happens during 
the initial transient, t <Tt ims . Here our simulations for 
P(t) reveal a quite complex behavior. After a very slow 
initial onset, a fast transient is found, which eventually 
ends in a local minimum marking the transition to the 
frozen relaxation regime, kth ~ 0. This transient regime 
does not extend into the nonadiabatic region. For in- 
stance, for T = 0.2hA/kB, it could only be observed for 
A/uj c > 0.2. 



V. CONCLUSIONS AND DISCUSSION 

Using the spin-boson model as a description for ET 
processes, we have calculated the thermal transfer rate 
k t h for two different reorganization energies, A = 10 A 
and A = 50 A, where A is twice the electronic coupling 
between the two redox states. Our real-time PIMC sim- 
ulations are able to cover the full range of temperatures 
and adiabaticity parameters A/u> c , where oj c denotes a 
typical prominent bath frequency. In this paper, we have 
confined ourselves to symmetric (unbiased) ET systems. 
The location of the calculated rates in parameter space 
are schematically shown in Fig. [l3|. The time-dependent 
function k(t), see Eq. (|9|), is a powerful means both to 
obtain k t h and to decide whether a rate description is ap- 
propriate in the first place. In the absence of electronic or 
vibrational coherence, this rate also describes the expo- 
nential decay of the electronic population P(t), provided 
the timescales for transient dynamics T trans and thermal 
relaxation r re i a x are sufficiently well separated. However, 
even in the absence of such a separation, P(t) can often 
be described by a rate with good accuracy. In this case 
we found that kth can be extracted from k(t) if the tran- 
sition between transient motion and the relaxation pro- 
cess can be clearly identified, this usually being the case 
if transient motion displays a pronounced peak. That 
there exist differences in the ET dynamics obtained from 
Pit) or k(t) should not come as a surprise, since these 
quantities correspond to different initial preparations. 

In the nonadiabatic regime, Eq. (|l2| ) captures the ther- 
mal transfer rate accurately for A/o; c <0.1. Even for 
A/uj c — 1 or 2, nonadiabatic theory typically overesti- 
mates k t h by only a small factor of the order 1 to 2, 
reflecting the increasing relevance of recrossing events ne- 
glected in Eq. (|l2|). Based on our simulations, we were 
also able to give precise bounds for the validity of the 
so-called scaling picture, see Eq. (|32|). The crossover to 
the adiabatic regime is then surprisingly well described 
by the famous Marcus theory, see Eq. ( |l4| ) . For high tem- 
peratures (within the classical regime), the solvent scale 
uj r in Eq. ( fl4"| ) is given by Eq. (fl6|). For lower tempera- 
tures, we identified an "extended" Marcus regime, where 
ui r is empirically given by Eq. ([m|). Eventually, for very 
low temperatures, T <hA/kB, ET dynamics is almost 
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completely restricted to transient motion, which cannot 
be explained by Marcus theory but requires a full dy- 
namical description. In this regime, analytical progress 
is most difficult, and our exact data reveal the presence 
of a "transient regime," where ET proceeds solely via 
this fast initial transient, while long-time relaxation is ex- 
tremely slow. Finally, we have shown that for A/lu c > 10, 
ET dynamics is fully coherent, and no rate description is 
possible. A comparison to the adiabatic description pl| 
reveals that the latter can only explain the dynamics of 
both P(t) and k(t) for times up to Tt ra ns- 

Our findings are summarized in the schematic "phase 
diagram" depicted in Fig. [IJ. Within large areas of pa- 
rameter space, ET can be described by an exponential 
decay after a transient timescale, t > r trans . Through- 
out this rate regime, the corresponding rate constant is 
very well captured by analytical approaches, namely the 
nonadiabatic approximation for small A/lu Ci and Mar- 
cus theory for intermediate-to-large A/cj c . However, a 
breakdown of the rate description can be observed in two 
different cases. In the strongly adiabatic regime, nuclear 
(vibrational) coherence causes oscillatory behavior. In 
addition, at very low temperatures within the crossover 
regime, relaxation is frozen in, and the transient motion 
dominates the ET process ( "transient regime" ) . 

Clearly, the computational method used in this paper 
has the potential to solve many other interesting prob- 
lems via real-time PIMC. In the context of ET dynamics, 
we have not reported results for the biased system, but 
plan to do so in the future. 



ACKNOWLEDGMENTS 

We wish to thank J. Ankerhold, M. Dikovsky, H. 
Grabert, A. Lucke, C.H. Mak, and J.T. Stockburger 
for valuable discussions. Financial support by the 
Volkswagen- Stiftung (No. 1/73 259) is acknowledged. 



APPENDIX A: COMPUTATIONAL DETAILS 

In this Appendix, some details and improvements 
over the approach of Ref. are summarized. Paths 
are blocked together in the following way. The 
discretized contour 7 is divided into L consecutive 
parts {zi,...,z qi },...,{z qL _ 1+ i,...,ZN}, corresponding 
to the different levels I = 1, . . . , L, dividing the spin paths 
{<7j} into L respective subsets {<Jj}i- Starting on the 
first level, for each configuration of the {<Jj}i>i the re- 
spective level- 1 block is formed by all possible configu- 
rations of the sub-path Switching to the second 
level, changing only the spins {cj}2 with all others fixed 
gives a level-2 block for each configuration of the {<Jj}i>2- 
Clearly, each level-2 block contains a whole set of level- 
1 blocks. This scheme is continued until on the highest 



level L, the final sampling of expectation values is carried 
out, with strongly reduced sign problem. 

To calculate the appropriate block averages, the total 
action must factorize with respect to the different levels. 
An appropriate partitioning of <&[{cr m , r]j, is given by 



^^£,jZj m + —cr' 1 (L2q+ m ,l + I<2qr+ro,2g+l) 



3=1 
91 <? 



k=2 j>k 



A 



'li 



^22 ' 



3=1 



^>i>i[{Vj,^j}i'>i} = X X& ( A i k & + iX jkVk) 

fe=9i_i+l j>k 

(Al) 



A22 \ ■* , 

j=qi-l+l 



(note Ajj = A 22 for 2 < j < q) with $[{cr m , r/j, ^}] = 
J2i=i ^i[{^m, rjji The contribution from the free 

propagator is split in the same manner. Here we in- 
clude the imaginary-time spins into the first level, since 
they contribute only weakly to the overall sign problem. 
Therefore they can be assigned to any level of real-time 
spins without significantly affecting the sign problem. 
Since the numerical effort to update spins increases for 
high levels, it saves computing time to put them into the 
lowest one. 

For a computation of C(t), a slight modification 
of previous MLB formulations is necessary since in 
Refs. |57]^60[| averages can be computed only on time 
slices belonging to the topmost level. This extension of 
MLB to a computation of an arbitrary n-point function 

A(h, . . . ,t n ) = (Aliaj}}) ~ Z- 1 J2 ^oMWj}} 

(A2) 

is described next. At first a dependence of A(t\, . . . , t n ) 
on all time slices seems to contradict the philosophy 
underlying MLB, namely tackling the sign problem by 
replacing the sampling process over all possible system 
paths by one over the level-L block only. However, de- 
composing A in the same way as $, 



A[W j }]=l[A l [{a j } 



v>i\ 



(A3) 
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block averages of Ai for I < L remove all depen- 
dence on lower-level spins. Therefore the remaining 
accumulation process can again be performed by sam- 
pling on the highest level only. These block aver- 
ages are calculated as modified "bonds" Bi [{vj}v>i\ — 
A i+i[{ a j}i'>i] B i[{ (J 3 }i'>i] where Bi denotes the "bond" 
of Ref. g, 

Bi[{v j }i<>t] = Cr_\ £ B l ^[{a j } v > l ]e- W ^M 

(A4) 

with the normalization C/_i = E{cr } ; "Pl—i [{ a j}i] of the 
weight for level-Z bonds, 

TVif-KM = \B l -i[{ t r j } h {* s }? +1 ,...,{* j }° L ] 

xe -m({^}i,{^}i+x>-,{^}l) . (A5) 

Sums run over all possible configuration of leveW spins 
{<Jj}i, and {cr°};/<z denotes their initial configuration at 
the start of the MC trajectory. Hence A(t\, . . . , t n ) can 
be written as 



influence functional 3>p[{<7j}'] is obtained from Eq. ( |30| ) 
by dropping all terms with imaginary-time spins a m and 
setting the t = spins <7i, a[ to +1. In addition 

C[fo}'] =iJ2^iQ(U ~ 1/2)A J ) - Q((J - 3/2)A,)} 
3=2 

(A8) 

includes a contribution describing the initial preparation 
of the bath To put P(t) into the form of Eq. (|A|), 

we rewrite Eq7(|26|) according to 

P(t) = U + x P -4^r) AM] , (A9) 



A(il,...,i n ) 
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W L [{a 3 }f 


)}| 



(A6) 



where the sums in the last expression run over K config- 
urations {cj }^ of the level- L spins distributed accord- 
ing to Vi\{&j\i}\- Here K is the usual MLB sampling 
number |57] , with typical values K k, 200 employed in 
the simulations. Therefore A(t\, . . . ,t n ) is accessible to 
a MLB scheme which not only soothes the sign problem 
due to interfering phase factors of different spin paths, 
but also takes care of alternating signs of the lower-level 
contributions to A. 

Finally we outline a scheme suitable for compensat- 
ing for different initial preparations, i.e. computing P(t) 
with the same MC trajectories as C(t). After discretizing 
the corresponding path integrals, P(t) can be written as, 
cp. Eq. pi), 



P[{°j}] 



with [note that E{ CTl)CT ; } E R „ } 1 = 2 2 • 2™] 

2 -m-2 z -l = Lm+2 p P [{a 3 Y] 



p[fo}] 



(A10) 



The effect of the different initial preparation is now 
absorbed by the factor pp/p. Consequently P(t) can 
now be computed as described above with = 
(T (j+iPp[{ 'j}']/p[{ 'j}]- The prefactor can be treated on 
the same footing. Thus, with the above extension in 
mind, P(t) can be calculated from the same MC tra- 
jectories used for sampling C{t) within our MLB-PIMC 
approach. 



P(t)=Z p 1 ]T a q+1 p P [{a 3 Y] 

p P [{a j }']=K(+l,<T2)K*(+l,a' 2 ) 
xe -*p[{^}']+C[{^}'] 



JY 



J\K{a l ,a l+ i) 



.i=2 



where the summation now only includes real-time spins 
<r(t > 0), i.e. {<7j}' = {(72, • • ■ ,&2q}- The corresponding 



[1] R.A. Marcus, J. Chem. Phys. 24, 966 (1956). 
[2] R.A. Marcus and N. Sutin, Biochim. Biophys. Acta 811, 
265 (1985). 

[3] A.M. Kuznetsov, Charge transfer in physics, chemistry 

and biology (Gordon and Breach, 1995). 
[4] H. Tributsch and L. Pohlmann, Science 279, 1891 (1998). 
[5] D. Chandler, in Liquids, Freezing and the Glass Transi- 
tion, ed. by D. Levesque, J. P. Hansen, and J. Zinn- Justin 
(Elsevier Science, North Holland, 1991), Les Houches 51, 
Part 1. 

[6] X. Song and A. A. Stuchebrukhov, J. Chem. Phys. 99, 
969 (1993). 

[7] U. Weiss, Quantum Dissipative Systems, Series in Mod- 
ern Condensed Matter Physics, Vol. 2 (World Scientific, 
Singapore, 1998). 

(A7) [8] R.P. Feynman and F.L. Vernon, Ann. Phys. (N.Y.) 24, 
118 (1963). 

[9] H. Grabert, P. Schramm, and G.-L. Ingold, Phys. Rep. 
168, 115 (1988). 



11 



[10] G.A. Voth, D. Chandler, and W.H. Miller, J. Chem. 

Phys. 93, 7009 (1989). 
[11] P. Hfinggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 

62, 251 (1990). 

[12] V.G. Levich, Adv. Electrochem. Electrochem. Eng. 4, 
249 (1965). 

[13] A. A. Ovchinnikov and M.Y. Ovchinnikova, Zh. Eksp. 
Teor. Fis. 56, 1278 (1969) [Sov. Phys. JETP 29, 688 
(1970)]. 

[14] MA. Vorotyntsev, R. Dogonadze, and A.M. Kuznetsov, 
Dokl. Akad. Nauk SSSR 195, 1135 (1970). 

[15] R. Egger, C.H. Mak, and U. Weiss, J. Chem. Phys. 100, 
2651 (1994). 

[16] B. Carmeli and D. Chandler, J. Chem. Phys. 82, 3400 

(1985); ibid. 89, 452 (1988). 
[17] J.N. Gehlen and D. Chandler, J. Chem. Phys. 97, 4958 

(1992). 

[18] A.I. Burshtein and Yu. Georgievskii, J. Chem. Phys. 100, 
7319 (1994). 

[19] L. Miihlbacher, A. Lucke, and R. Egger, J. Chem. Phys. 

110, 5851 (1999). 

[20] M.J. Hornbach and Yu. Dakhnovskii, J. Chem. Phys. 

111, 5073 (1999). 

[21] A. Garg, J.N. Onuchic, and V. Ambegaokar, J. Chem. 

Phys. 83, 4491 (1985). 
[22] F. Lesage and H. Saleur, Phys. Rev. Lett. 80, 4370 

(1998). 

[23] A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P.A. 
Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 
1 (1987). 

[24] R. Egger, C.H. Mak, and U. Weiss, Phys. Rev. E 50, 655 

(1994) . 

[25] M. Grifoni, E. Paladino, and U. Weiss, Eur. Phys. J. B 
10, 719 (1999). 

[26] R. Silbey and R.A. Harris, J. Chem. Phys. 80, 2615 
(1984). 

[27] W.Th. Pollard, A.K. Felts, and R.A. Friesner, Adv. 

Chem. Phys. 93, 77 (1996). 
[28] G. Stock, Phys. Rev. E 51, 3038 (1995). 
[29] G. Stock, J. Chem. Phys. 103, 1561 (1995). 
[30] G. Stock and M. Thoss, Phys. Rev. Lett. 78, 578 (1997). 
[31] X. Sun, H. Wang, and W.H. Miller, J. Chem. Phys. 109, 

7064 (1998). 

[32] H. Wang, X. Song, D. Chandler, and W.H. Miller, J. 

Chem. Phys. 110, 4828 (1999). 
[33] A. Lucke, J. Stockburger, and C.H. Mak, J. Chem. Phys. 

Ill, 10843 (1999). 
[34] A. A. Golosov, R.A. Friesner, and P. Pechukas, J. Chem. 

Phys. 110, 138 (1999). 
[35] A. A. Golosov, R.A. Friesner, and P. Pechukas, J. Chem. 

Phys. 112, 2095 (2000). 
[36] J. Casado-Pascual, C. Denk, M. Morillo, and R.I. Cukier, 

J. Chem. Phys. 113, 11176 (2000). 
[37] H. Wang, M. Thoss, and W.H. Miller, J. Chem. Phys. 

115, 2979 (2001); M. Thoss, H. Wang, and W.H. Miller, 

J. Chem. Phys. 115, 2991 (2001). 
[38] H. Wang, J. Chem. Phys. 113, 9948 (2000). 
[39] N. Makri and D.E. Makarov, J. Chem. Phys. 102, 4600 

(1995) . 

[40] M. Winterstetter and W. Domcke, Chem. Phys. Lett. 
236, 445 (1995). 



[41] J. Stockburger and C.H. Mak, Phys. Rev. Lett. 80, 2657 
(1998). 

[42] J. Stockburger and H. Grabert, Phys. Rev. Lett. 88, 
170407 (2002). 

[43] T.A. Costi and C. Kieffer, Phys. Rev. Lett. 76, 1683 
(1996). 

[44] T.A. Costi, Phys. Rev. Lett. 80, 1038 (1998). 
[45] M. Keil and H. Schoeller, Phys. Rev. B 63, 180302 
(2001). 

[46] Quantum Monte Carlo Methods in Condensed Matter 

Physics, ed. by M. Suzuki (World Scientific, Singapore, 

1993), and references therein. 
[47] See, for instance, E.Y. Loh Jr., J. Gubernatis, R.T. 

Scalettar, S.R. White, D.J. Scalapino, and R.L. Sugar, 

Phys. Rev. B 41, 9301 (1990). 
[48] S. Chakravarty and J. Rudnick, Phys. Rev. Lett. 75, 501 

(1995). 

[49] K. Volker, Phys. Rev. B 58, 1862 (1998). 

[50] D. Bailey, M. Hurley, and H.K. McDowell, J. Chem. 

Phys. 109, 8262 (1998). 
[51] E.C. Behrman and P.G. Wolynes, J. Chem. Phys. 83, 

5863 (1985). 

[52] C.H. Mak and D. Chandler, Phys. Rev. A 44, 2352 
(1991). 

[53] R. Egger and C.H. Mak, Phys. Rev. B 50, 15210 (1994). 

[54] V.S. Filinov, Nucl. Phys. B 271, 717 (1986). 

[55] J.D. Doll, T.L. Beck, and D.L. Freeman, J. Chem. Phys. 

89, 5753 (1988). 
[56] C.H. Mak and R. Egger, Adv. Chem. Phys 93, 39 (1996), 

and references therein. 
[57] C.H. Mak, R. Egger, and H. Weber-Gottschick, Phys. 

Rev. Lett. 81, 4533 (1998). 
[58] C.H. Mak and R. Egger, J. Chem. Phys. 110, 12 (1999). 
[59] R. Egger, L. Miihlbacher, and C.H. Mak, Phys. Rev. E 

61, 5961 (2000). 
[60] M. Dikovsky and C.H. Mak, Phys. Rev. B 63, 235105 

(2001). 

[61] We consider only a > 1/2. ET dynamics is normally 
incoherent with respect to electronic motion, while for 
a < 1/2, spin-boson dynamics shows electronic coherence 

§■ 



12 



Im: 








1 .5 




j = 2q+l 


i 


. ' Re: 

j = q+l 


o 




j = 2q + r+ 1 






^0.5 




g = -mfi 









FIG. 1. KadanofF-Baym contour 7 in the complex-time 
plane and a possible discretization. 
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FIG. 2. Thermal transfer rate k t h for a — 1/2, uj c = 50A 
(diamonds, left axis) and uj c — 10A (squares, right axis). 
Filled symbols refer to values obtained from k(t), open ones 
to those obtained from P(t). Vertical lines indicate error bars 
from statistical sampling uncertainties. The error in fitting 
P(t) to an exponential decay is of the order of the symbol 
size. The dashed-dotted curves represent the nonadiabatic 
result (M). 
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FIG. 3. Same as Fig. Q but as a function of w c /A for fixed 
T — hA/ks- The dotted curve denotes the rate according to 
the scaling limit, while the dashed-dotted one gives Eq. (EOh . 
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FIG. 4. Thermal transfer rate kth f° r a — 1 with A = 50A 
(diamonds, left axis) and A = 10 A (squares, right axis). The 
curves give Eq. (^l]) . 




FIG. 5. P(t) and k(t) for lu c /A = 25, A = 50A, 
and T = TiA/kB- The dashed-dotted line represents 
exp(— 0.0042A t) — 0.008, in agreement with the plateau value 
of k(t). 
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FIG. 6. Same as Fig. | but for A = 10A, A/uj c = 1, 
and T = 3.333ftA/fc_B- The dashed-dotted line represents 
exp(-0.057A t) - 0.12. 
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FIG. 7. ET rate kth for A = 10A, with lj c = A (diamonds) 
and oJc = 0.5 A (circles). The dashed and solid curves repre- 
sent the corresponding Marcus rates ( p^ ) with u r = lu c /2. 
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FIG. 8. Same as Fig. ^| but for lower temperatures. Data 
for T = O.lhA/kB (not shown) are virtually identical to those 
for T = 0.2hA/kB , indicating that the zero-temperature limit 
has been reached. 
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FIG. 9. Same as Fig. g but for A/w c = 2, T = 2hA/k B 
(circles), T = hA/k B (diamonds), T = 0.6hA/k B (squares) 
and T = 0.2hA/k B (triangles). The left side of the graph 
shows the initial transient dynamics. The box in the top left 
corner of the upper right graph represents the correspond- 
ing scale of the upper left graph. The dotted [dashed, long 
dashed, dot-dashed] curve represents the adiabatic approxi- 
mation for T = 2hA/k B [hA/k B , 0.6hA/k B , 0.2hA/k B }. 
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FIG. 10. P(t) and k(t) for T = hA/k B , A/lu c = 5 (left) 
and A/uj c — 10 (right). Dot-dashed curves represent the adi- 
abatic prediction. 
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FIG. 11. Same as Fig. [lo[ but for T = 10hA/k B and 
(from left to right) A/lo c = 5, 10 and 20. 
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FIG. 12. Thermal transfer rate kth as a function of u c 
for T = WhA/k B (circles), T = 2hA/k B (triangles) and 
T — HA/kB (diamonds). The dotted [dashed, dashed-dotted] 
lines refer to the Marcus rate with uj r according to Eq. (^) 
and q(T = WhA/k B ) = 0.0002 [q(T = 2hA/k B ) = 0.21, 
q(T = TiA/ks) = 0.26]. The exact rates exceed the nonadi- 
abatic (maximum) Marcus rate for u c >, 5A [2A, 1.5A]. The 
inset shows kth at T = 2ftA/fcs and u c = 0.5 A as a func- 
tion of the reorganization energy A, where the dashed curve 
corresponds to the main graph. 
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FIG. 13. Location of parameter values where simula- 
tions have been carried out. The reorganization energies are 
A = 50A for A/u c = 0.02, 0.04, and A = 10A otherwise. Re- 
sults following the nonadiabatic prediction ( |l2| ) are denoted 
by circles, those following the (classical or extended) Marcus 
formula ([HJ) by squares. Transient behavior is marked by 
triangles, and oscillatory (adiabatic) dynamics by diamonds. 
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FIG. 14. Schematic ET dynamics "phase diagram." To 
obtain quantitative numbers for the crossover lines, refer to 
Fig. §. 



